home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / basic / _rational.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  7.4 KB  |  356 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  _rational.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16. #include <LEDA/basic.h>
  17. #include <LEDA/rational.h>
  18. #include <math.h>
  19. #include <ctype.h>
  20.  
  21.  
  22.  
  23. // hidden functions
  24.  
  25. LedaRational& LedaRational::normalize()
  26. // divide numerator and denominator by their greatest common divisor
  27.   { // den is assumed to be nonzero and positive
  28.     if (num == den) { num = den = 1; return (*this); }
  29.     if (-num == den) { num = -1; den = 1; return (*this); }
  30.     Int ggt = bgcd(num, den);
  31.     if (!Ieq1(ggt)) {
  32.       num /= ggt;
  33.       den /= ggt;
  34.     };
  35.     return (*this);
  36.   };
  37.  
  38.  
  39.  
  40. // Constructors
  41.  
  42. LedaRational::LedaRational(double x)
  43. // from the GNU-C++ library (MACHINE-DEPENDENT)
  44.   { num = 0; den = 1;
  45.  
  46.     if (x != 0.0)
  47.     { int neg = (x < 0);
  48.       if (neg) x = -x;
  49.  
  50.       const unsigned shift = 15;     // a safe shift per step
  51.       const double width = 32768;    // = 2^shift
  52.       const int maxiter = 20;        // ought not be necessary, but just in case,
  53.                                      // max 300 bits of precision
  54.       int expt;
  55.       double mantissa = frexp(x, &expt);
  56.       long exponent = expt;
  57.       double intpart;
  58.       int k = 0;
  59.       while (mantissa != 0.0 && k++ < maxiter) {
  60.         mantissa *= width; // shift double mantissa
  61.         mantissa = modf(mantissa, &intpart);
  62.         num <<= shift;
  63.         num += (long)intpart;
  64.         exponent -= shift;
  65.       }
  66.       if (exponent > 0)
  67.         num <<= (unsigned)exponent;
  68.       else if (exponent < 0)
  69.         den <<= (unsigned)(-exponent);
  70.       if (neg)
  71.         num.negate();
  72.     } // if (x != 0) then
  73.     (*this).normalize();
  74.   };
  75.  
  76.  
  77. LedaRational::LedaRational(int n, int d)
  78.   { if (d == 0) {
  79.       error_handler(1,"Zero denominator!");
  80.     }
  81.     else {
  82.       num = Int(n);
  83.       den = Int(d);
  84.       if (d < 0) {
  85.         num.negate();
  86.         den.negate();
  87.       }
  88.       (*this).normalize();
  89.     }
  90.   };
  91.  
  92.  
  93. LedaRational::LedaRational(const Int& n, const Int& d)
  94.   { if (Ieq0(d)) {
  95.       // d == 0
  96.       error_handler(1,"Zero denominator!");
  97.     }
  98.     else {
  99.       num = n;
  100.       den = d;
  101.       if (Ilt0(d)) {
  102.         // d < 0
  103.         num.negate();
  104.         den.negate();
  105.       }
  106.       (*this).normalize();
  107.     }
  108.   };
  109.  
  110.  
  111.  
  112. // Arithmetic Operators
  113.  
  114. LedaRational& LedaRational::operator+= (const LedaRational& r)
  115.   { num = num * r.den + r.num * den;
  116.     den *= r.den;
  117.     return (*this).normalize();
  118.   };
  119.  
  120. LedaRational& LedaRational::operator-= (const LedaRational& r)
  121.   { num = num * r.den - r.num * den;
  122.     den *= r.den;
  123.     return (*this).normalize();
  124.   };
  125.  
  126. LedaRational& LedaRational::operator*= (const LedaRational& r)
  127.   { num *= r.num;
  128.     den *= r.den;
  129.     return (*this).normalize();
  130.   };
  131.  
  132. LedaRational& LedaRational::operator/= (const LedaRational& r)
  133.   { if (Ieq0(r.num)) {
  134.       // r == 0
  135.       error_handler(1,"Division by 0!");
  136.     }
  137.     else {
  138.       // r.num != 0
  139.       num *= r.den;
  140.       den *= r.num;
  141.       if (Ilt0(den)) {
  142.         num.negate();
  143.         den.negate();
  144.       }
  145.     }
  146.     return (*this).normalize();
  147.   };
  148.  
  149.  
  150.  
  151. // Assignment Operator
  152.  
  153. LedaRational& LedaRational::operator= (const LedaRational& r)
  154.   { if (this == &r) return *this; // to handle r = r correctly
  155.     num = r.num;
  156.     den = r.den;
  157.     return *this;
  158.   };
  159.  
  160.  
  161.  
  162.  
  163. // some useful member-functions
  164.  
  165. void LedaRational::invert()
  166.   { if (Ieq0(num)) {
  167.       error_handler(1,"Zero denominator!");
  168.     }
  169.     else {
  170.       Int tmp = num;
  171.       num = den;
  172.       den = tmp;
  173.       if (Ilt0(den)) {
  174.         num.negate();
  175.         den.negate();
  176.       }
  177.     }
  178.   };
  179.  
  180. LedaRational LedaRational::inverse()
  181.   { if (Ieq0(num)) {
  182.       error_handler(1,"Zero denominator!");
  183.       return (*this);
  184.     }
  185.     else {
  186.       LedaRational tmp(den,num);
  187.       if (Ilt0(num)) {
  188.         (tmp.num).negate();
  189.         (tmp.den).negate();
  190.       }
  191.       return tmp;
  192.     }
  193.   };
  194.  
  195.  
  196.  
  197. // friend functions
  198.  
  199. int compare(const LedaRational& x, int y)
  200.   { int
  201.       xsign = sign(x.num),
  202.       ysign;
  203.     if (y == 0) ysign = 0;
  204.     else ysign = (y > 0) ? 1 : -1;
  205.     if (xsign == 0) { return -ysign; }
  206.     if (ysign == 0) { return xsign; }
  207.     // now (x != 0) && (y != 0)
  208.     int
  209.       diff = xsign - ysign;
  210.     if (diff == 0) {
  211.       Int
  212.         leftop = x.num,
  213.         rightop = Int(y) * x.den;
  214.       if (leftop < rightop) { return -1; }
  215.       else { return (leftop > rightop); }
  216.     }
  217.     else return diff;
  218.   };
  219.  
  220.  
  221. int compare(int x, const LedaRational& y)
  222.   { int
  223.       ysign = sign(y.num),
  224.       xsign;
  225.     if (x == 0) xsign = 0;
  226.     else xsign = (x > 0) ? 1 : -1;
  227.     if (xsign == 0) { return -ysign; }
  228.     if (ysign == 0) { return xsign; }
  229.     // now (x != 0) && (y != 0)
  230.     int
  231.       diff = xsign - ysign;
  232.     if (diff == 0) {
  233.       Int
  234.         leftop = Int(x) * y.den,
  235.         rightop = y.num;
  236.       if (leftop < rightop) { return -1; }
  237.       else { return (leftop > rightop); }
  238.     }
  239.     else return diff;
  240.   };
  241.  
  242.  
  243. // other useful friend functions
  244.  
  245. LedaRational pow(const LedaRational& r, int l)
  246. // no need to normalize since num and den are relatively prime
  247.   { LedaRational mul(r), result(1,1);
  248.     if (l < 0) {
  249.       mul = mul.inverse();
  250.       l = -l;
  251.     }
  252.     for (int i = 1; i <= l; i++) {
  253.       result.num *= mul.num;
  254.       result.den *= mul.den;
  255.     }
  256.     return result;
  257.   };
  258.  
  259. LedaRational pow(const LedaRational& r, Int I)
  260. // no need to normalize since num and den are relatively prime
  261.   { LedaRational mul(r), result(1,1);
  262.     if (Ilt0(I)) {
  263.       mul = mul.inverse();
  264.       I.negate();
  265.     }
  266.     for (Int i = 1; i < I; i++) {
  267.       result.num *= mul.num;
  268.       result.den *= mul.den;
  269.     }
  270.     return result;
  271.   };
  272.  
  273. LedaRational::operator double() const
  274.   { Int
  275.       numvar = num,
  276.       denvar = den;
  277.  
  278.     if (numvar == Int(0)) { return 0; }
  279.     
  280.     const Int MDP = 1000000;    // my_double_precision
  281.     long s = 0;
  282.     Int quot = (numvar / denvar); // integer quotient
  283.  
  284.     while (abs(quot) < MDP) {
  285.       numvar *= 10;
  286.       s++;
  287.       quot = (numvar / denvar);
  288.     }
  289.     // |quot| > MDP
  290.     // num_new == 10^s * num_old
  291.  
  292.     while (abs(quot) > MDP) {
  293.       quot /= 10;
  294.       s--;
  295.     }
  296.     // MDP/10 < |quot| < MDP
  297.     // num_old/den_old == quot * 10^{-s}
  298.  
  299.     double result = (double)longasI(quot);
  300.     // transform Int into double via long
  301.  
  302.     if (s >= 0) {
  303.       for (int i = 0; i < s; i++) { result /= 10; };
  304.     }
  305.     else {
  306.       for (int i = 0; i > s; i--) { result *= 10; };
  307.     }
  308.     return result;
  309.   }; 
  310.  
  311. Int floor(const LedaRational& r)
  312.   { Int x, y;
  313.     Idiv (x, y, r.num, r.den);
  314.     if ((Ilt0(r.num)) && (!Ieq0(y))) x--;
  315.     return x;
  316.   };
  317.  
  318. Int ceil(const LedaRational& r)
  319.   { Int x, y;
  320.     Idiv (x, y, r.num, r.den);
  321.     if ((Ige0(r.num)) && (!Ieq0(y))) x++;
  322.     return x;
  323.   };
  324.  
  325. Int round(const LedaRational& r)
  326.   { Int rem, quot;
  327.     Idiv(quot, rem, r.num, r.den);
  328.     rem <<= 1;
  329.     if (rem >= r.den) {
  330.       if (sign(r.num) >= 0) { quot++; }
  331.       else { quot--; }
  332.     }
  333.     return quot;
  334.   }
  335.  
  336. istream& operator>> (istream& in, LedaRational& r)
  337.   { // Format: "r.num / r.den"
  338.     Int rx, ry;
  339.     char c;
  340.     do in.get(c); while (isspace(c));
  341.     in.putback(c);
  342.  
  343.     in >> rx;
  344.    
  345.     do in.get(c); while (isspace(c));
  346.     if (c != '/') { error_handler(1,"/ expected"); }
  347.  
  348.     do in.get(c); while (isspace(c));
  349.     in.putback(c);
  350.    
  351.     in >> ry;
  352.     r = LedaRational(rx,ry);
  353.     // to guarantee the value is normalized, denominator is nonzero ...
  354.     return in;
  355.   };
  356.